From raw data to start of analysis

The data sources as well as some preliminary data cleaning steps are inspired by a FiveThirtyEight article from 2014: The Economic Guide To Picking A College Major.

Our primary data source is the American Community Survey (ACS) Public Use Microdata Sample (PUMS) conducted by the United States Census Bureau. The ACS conducts surveys of around 3.5 million households each year.

We are using the multiyear PUMS file covering census data for the years 2012 to 2016, which contains records on ~16 million individual residents. It can be downloaded as zip file by selecting “United States Population Records” from this website. The zip file contains four large csv files which together cover the complete dataset. Additional information on the PUMS files is provided on the US Census Bureau’s website: they provide a general overview (link to PDF) as well as a data dictionary (link to PDF).

In the first data cleaning step, we remove individuals who did not graduate college (see R script here) and write the remaining 3,489,070 entries to a file called us16_filtered.csv.

Next up, we add information about college majors by merging our dataset with a lookup file called majors-list.csv prepared by FiveThirtyEight (GitHub link). This lookup file matches the FOD1P codes used in the PUMS to actual names of the corresponding college majors. Additionally, these college majors are grouped into categories as used in a major 2011 study 44 Carnevale et al., “What’s It Worth?: The Economic Value of College Majors.” Georgetown University Center on Education and the Workforce, 2011. Website.

You can find the R script we used to merge datasets on GitHub. This script also converts all wages reported before 2016 to inflation-adjusted wages equivalent to 2016 dollars, before writing selected columns for the 3,489,070 entries to a file called selected_vars.csv.

Intermission: what is it we actually want to do?

Before moving on with the data preparation steps, it is probably a good idea to set the stage regarding what we want to model and why. Hopefully, this will make some of the following sections easier to follow.

As touched upon in the introduction, our intention here is to run relatively straightforward hierarchical regression models in order to model wages across different ages. Doing this in a Bayesian setting of course allows us to draw posterior predictive samples from the fitted models and we can use these posterior draws to answer some questions about the data.

For instance, you will see that the wage gap in absolute dollar values is substantially smaller when looking at education graduates as opposed to business graduates. However, wages are generally lower in the education group, so the observation about the wage gap might just be a side effect of this general difference in wages. Drawing posterior predicitve samples allows us to check if the difference in the wage gap holds up when looking at relative values rather than absolute ones.

Selecting the data subset

The file selected_vars.csv is the starting point for our actual analysis:

library(readr)
selected_vars <- read_csv("data/selected_vars.csv")
selected_vars

## # A tibble: 3,489,070 x 9
##     SERIALNO    ST  AGEP   SEX FOD1P Major  Major_Category real_total_earn
##        <dbl> <int> <int> <int> <int> <chr>  <chr>                    <dbl>
##  1   2.01e¹²    13    69     1  4101 PHYSI… Industrial Ar…           36961
##  2   2.01e¹²    13    53     2  2304 ELEME… Education                41185
##  3   2.01e¹²    13    34     2  1100 GENER… Agriculture &…            4435
##  4   2.01e¹²    13    56     1  6203 BUSIN… Business                 13728
##  5   2.01e¹²     6    49     1  2414 MECHA… Engineering              72866
##  6   2.01e¹²     6    48     2  3600 BIOLO… Biology & Lif…          211206
##  7   2.01e¹²     6    50     2  2300 GENER… Education                84482
##  8   2.01e¹²     5    67     2  6004 COMME… Arts                     10455
##  9   2.01e¹²     6    59     1  5003 CHEMI… Physical Scie…               0
## 10   2.01e¹²     6    54     2  2901 FAMIL… Industrial Ar…               0
## # ... with 3,489,060 more rows, and 1 more variable: real_wage <dbl>

It contains nine variables:

• SERIALNO: ID number for each entry
• ST: two-digit code for US states
• AGEP: age of survey respondent (at time of answering)
• SEX: survey respondent’s gender
• FOD1P: four-digit code for college major
• Major: corresponding college major (173 majors in total)
• Major_Category: college major category (16 categories in total)
• real_total_earn: inflation-adjusted reported yearly total earnings (in 2016 USD)
• real_wage: inflation-adjusted reported yearly wages (in 2016 USD)

While it would of course be most appropriate to conduct an analysis on the complete dataset, we chose to focus only on a select subset of entries in order to keep both computation times and model interpretation manageable. Our chosen subset consists of entries from two states and three major categories.

We choose New York and Texas as two similarly large states which should serve as decent stand-ins for politically opposed blue states (NY) and red states (TX). The major categories we focus on are business, education and health. We select the data subset with the following R code:

library(dplyr)
nytx_df <- selected_vars %>% 
  filter(ST %in% c(36, 48),
         Major_Category %in% c("Education", "Health", "Business"),
         real_wage > 0) %>% 
  mutate(state = if_else(ST == 36, "NY", "TX"),
         sex = if_else(SEX == 1, "male", "female"))

Note that in the above code chunk, we also remove respondents who reported no wage and recode ST and SEX into the more explicit variables state and sex.

The resulting subset contains 148,490 entries, which can be broken down as follows:

library(tidyr)
nytx_df %>% 
  group_by(Major_Category, state, sex) %>% 
  summarise(entries = n()) %>% 
  spread(key = "sex", value = "entries") %>% 
  mutate(ratio = female / male,
         total = male + female) %>% 
  knitr::kable(., digits = 2, format.args = list(big.mark = ","),
               col.names = c("major category", "state", "female", "male",
                             "ratio (f/m)", "total"))

major category	state	female	male	ratio (f/m)	total
Business	NY	15,789	19,265	0.82	35,054
Business	TX	18,888	24,829	0.76	43,717
Education	NY	15,422	4,996	3.09	20,418
Education	TX	17,681	5,024	3.52	22,705
Health	NY	11,075	2,507	4.42	13,582
Health	TX	10,417	2,597	4.01	13,014

We can see that both total number and gender ratio of entries for the two states are roughly comparable across the three major categories.

Data summaries

In order to get an overview of the general trend in wages as a function of age, we plot the mean and variation of yearly wages across the age range from 20 to 70:

library(ggplot2)
theme_set(custom_ggtheme)

nytx_df %>% 
  filter(AGEP >= 20,
         AGEP <= 70) %>% 
  group_by(AGEP, sex, Major_Category, state) %>% 
  summarise(mean_wage = mean(real_wage, na.rm = T),
            sd_wage = sd(real_wage, na.rm = T)) %>% 
  ggplot(., aes(x = AGEP)) +
  geom_line(aes(y = mean_wage, colour = sex)) +
  geom_line(aes(y = mean_wage + 2 * sd_wage, colour = sex), linetype = "dashed") +
  ggtitle("Data summary: yearly wages (mean + 2 standard deviations)") +
  xlab("age") +
  ylab("yearly wages in USD") +
  facet_grid(state ~ Major_Category)

We can see differences in wages, slopes as well as in the gender wage gap when comparing the different panels. Nothing here is too surprising: business graduates have the highest earning potential, followed by graduates in the field of health, and education graduates coming in last. Furthermore, it appears that wages grow steadily through people’s thirties and forties until growth starts to stall around the age of fifty. This observation is fairly consistent across all panels.

Looking at the gender wage gap, the gap is consistently smallest within education graduates (at least in terms of absolute numbers). While the gap within health and business ends up being roughly similar in older age cohorts, the way there is different: women with a business background seem to earn less than men right out of the gate, whereas wages only start to separate around age 35 for health graduates.

On a more technical note, we can also see that there is a lot of variation in the wages within each age cohort. This is not unexpected, however it has implications for our model choice: If we were to go ahead and use a regression model that directly models wage (with Gaussian noise) as a function of age, it is very likely that we would end up with negative wages when drawing posterior predictive samples from that model.

In order to avoid this, we log-transform the wage numbers before running our models. This will ensure that our posterior predictive samples are positive (after exponentiating them) – and has the convenient side-effect of lowering the impact of high outlier wages.

model_df1 <- nytx_df %>% 
  mutate(log_wage = log(real_wage))

While using Gaussian distributions to describe the wage distributions at each age point is of course still an approximation, it is not a terribly outrageous one after log-transforming the wages. We can see this in the violin plots of select age cohorts (we only show eight because it would start to get cramped with more age cohorts):

model_df1 %>% 
  filter(AGEP %in% seq(25, 60, 5)) %>% 
  ggplot(., aes(x = factor(AGEP), y = log_wage)) +
  geom_violin() +
  ggtitle("Data summary: log-transformed yearly wages in eight age cohorts") +
  xlab("age") +
  ylab("log-transformed yearly wages") +
  facet_grid(state ~ Major_Category)

lin_df <- model_df1 %>% 
  filter(AGEP >= 25,
         AGEP <= 50)

set.seed(4242)
prototype_lin_df <- lin_df %>% 
  sample_n(size = 5000)

Model fitting & convergence diagnostics

We can now go ahead and fit a hierarchical generalised linear model using the stan_glmer() function from rstanarm. Log-transformed wages will be fitted as a linear function of age with Gaussian noise. We want to allow for the slope of this function to be different between the genders, so we add sex as a grouping term. In addition, we specify the hierarchical groupings Major_Category and state.

Sampling takes place with rstanarm’s default settings: four chains with 4,000 samples, 2,000 of which are warmup samples. Default priors are N(0, 2.5) for regression coefficients, N(0, 10) for intercepts and exponential(1) for sigma parameters of the Gaussian noise. Note that the prior parameters apply after the regression’s predictors are centered.

library(rstan)
library(rstanarm)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores() - 1)

lin_fit <- stan_glmer(log_wage ~ (AGEP | sex) + 
                        (1 | Major_Category) + (1 | state),
                      data = prototype_lin_df, family = gaussian(), 
                      adapt_delta = 0.99)

Let’s have a look at the model fit summary, including parameter estimates and convergence diagnostics:

knitr::kable(lin_fit$stan_summary, digits = 2)

	mean	se_mean	sd	2.5%	10%	25%	50%	75%	90%	97.5%	n_eff	Rhat
(Intercept)	10.03	0.02	0.58	8.76	9.49	9.81	10.04	10.26	10.53	11.06	599.10	1.01
b[(Intercept) Major_Category:Business]	0.09	0.01	0.30	-0.44	-0.13	-0.02	0.08	0.17	0.30	0.67	967.76	1.00
b[(Intercept) Major_Category:Education]	-0.17	0.01	0.30	-0.67	-0.40	-0.28	-0.18	-0.08	0.05	0.40	972.32	1.00
b[(Intercept) Major_Category:Health]	0.11	0.01	0.30	-0.39	-0.11	0.00	0.10	0.20	0.33	0.70	951.42	1.00
b[(Intercept) Major_Category:_NEW_Major_Category]	0.00	0.01	0.49	-0.91	-0.38	-0.16	0.00	0.15	0.38	0.90	2876.48	1.00
b[(Intercept) state:NY]	0.06	0.02	0.45	-0.76	-0.23	-0.04	0.05	0.16	0.40	0.92	554.35	1.01
b[(Intercept) state:TX]	-0.05	0.02	0.45	-0.87	-0.34	-0.15	-0.05	0.05	0.29	0.82	554.60	1.01
b[(Intercept) state:_NEW_state]	0.00	0.01	0.67	-1.33	-0.45	-0.14	0.00	0.14	0.47	1.30	3276.54	1.00
b[(Intercept) sex:female]	0.07	0.01	0.19	-0.21	-0.06	-0.01	0.04	0.14	0.26	0.51	1184.33	1.00
b[AGEP sex:female]	0.02	0.00	0.00	0.01	0.01	0.01	0.02	0.02	0.02	0.02	3118.43	1.00
b[(Intercept) sex:male]	-0.09	0.01	0.19	-0.51	-0.29	-0.15	-0.05	0.00	0.05	0.22	1351.24	1.00
b[AGEP sex:male]	0.03	0.00	0.00	0.03	0.03	0.03	0.03	0.03	0.03	0.04	2581.58	1.00
b[(Intercept) sex:_NEW_sex]	0.01	0.00	0.29	-0.52	-0.21	-0.07	0.00	0.08	0.24	0.59	3312.71	1.00
b[AGEP sex:_NEW_sex]	0.00	0.00	0.11	-0.24	-0.09	-0.03	0.00	0.03	0.08	0.22	2025.84	1.00
sigma	0.93	0.00	0.01	0.91	0.92	0.93	0.93	0.94	0.94	0.95	4000.00	1.00
Sigma[Major_Category:(Intercept),(Intercept)]	0.23	0.02	0.93	0.01	0.01	0.03	0.06	0.16	0.47	1.60	1464.14	1.00
Sigma[state:(Intercept),(Intercept)]	0.44	0.05	1.48	0.00	0.00	0.01	0.05	0.26	1.00	3.44	1017.29	1.00
Sigma[sex:(Intercept),(Intercept)]	0.09	0.01	0.35	0.00	0.00	0.00	0.02	0.07	0.20	0.64	2120.05	1.00
Sigma[sex:AGEP,(Intercept)]	0.00	0.00	0.03	-0.05	-0.01	0.00	0.00	0.00	0.01	0.05	2242.27	1.00
Sigma[sex:AGEP,AGEP]	0.01	0.00	0.03	0.00	0.00	0.00	0.00	0.01	0.03	0.10	1192.68	1.00
mean_PPD	10.87	0.00	0.02	10.83	10.85	10.86	10.87	10.88	10.89	10.90	4000.00	1.00
log-posterior	-6775.87	0.12	3.93	-6784.42	-6781.16	-6778.33	-6775.53	-6773.14	-6771.06	-6768.95	1088.16	1.00

All Rhat values are below 1.1, giving an indication that the chains have converged. We can also asses convergence by inspecting the trace plots of the sampling chains:

library(bayesplot)
theme_set(custom_ggtheme)

mcmc_trace(as.array(lin_fit), np = nuts_params(lin_fit))

DIVERGENCES? (Tried to increase delta more and more. no dice so far)

Predictive performance: PSIS-LOO

We can use Pareto smoothed importance sampling (PSIS) leave-one-out (LOO) cross-validation (CV) to assess the predictive performance of the linear model. This is done by applying LOO-CV to the log-likelihood matrix of the fitted model. We can then plot the estimated shape parameters k of the generalised Pareto distribution in order to check the reliability of the LOO expected log pointwise predictive (elpd) estimate:

loo_lin <- loo(log_lik(lin_fit))
par(bg = '#fffff8')
plot(loo_lin)

All k values are below 0.5, indicating quick convergence and solid predicite performance 55 Vehtari, Aki, Andrew Gelman, and Jonah Gabry. “Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.” Statistics and Computing 27.5 (2017): 1413-1432..

Posterior predictions

With the linear model fitted, we can now draw posterior predictive samples to answer some of our questions around wages and the gender wage gap.

First, let’s have a look at the predicted wages. We show the median and central 95% interval for wages in each age cohort:

post_sample_grid <- expand.grid(AGEP = 25:50, 
                                sex = c("male", "female"), 
                                Major_Category = c("Education", "Health", "Business"), 
                                state = c("NY", "TX"),
                                KEEP.OUT.ATTRS = F,
                                stringsAsFactors = F)

log_post_lin <- posterior_predict(lin_fit, post_sample_grid, seed = 42)
exp_post_lin <- exp(log_post_lin)

post_df <- bind_cols(post_sample_grid, as.data.frame(t(exp_post_lin)))

post_long <- post_df %>% 
  gather(key = "sample_no", value = "predicted_earning", V1:V4000)

ci_post_pred <- post_long %>% 
  group_by(AGEP, sex, Major_Category, state) %>% 
  summarise(q_lower = quantile(predicted_earning, 0.025, names = F),
            median = quantile(predicted_earning, 0.5, names = F),
            q_upper = quantile(predicted_earning, 0.975, names = F),
            q90 = quantile(predicted_earning, 0.9, names = F)) %>% 
  mutate(model = "linear")

ci_post_pred %>% 
  ggplot(., aes(x = AGEP)) +
  geom_line(aes(y = median, colour = sex)) +
  geom_line(aes(y = q_lower, colour = sex), linetype = "dashed") +
  geom_line(aes(y = q_upper, colour = sex), linetype = "dashed") +
  ggtitle("Posterior prediction (linear model): yearly wages (median & central 95% interval)") +
  xlab("age") +
  ylab("yearly wages in USD") +
  facet_grid(state ~ Major_Category)

At a first glance, we can already see that the model is not entirely unreasonable: it has generally replicated the trends we have seen in the real data, including the existence of a gender wage gap.

In order to show the context of real data more explicitly, we can plot the predictions together with the real data. Note that his time around, we plot the data median and the 90% quantile at each age point (instead of mean & two standard deviations earlier). Due to high wages being masked in our original dataset, the real data’s 97.5% wage quantiles look pretty peculiar (they form a horizontal line across all ages). As a consequence, the real data may look different in this plot compared to the earlier one.

data_intervals_median <- model_df1 %>% 
  filter(AGEP >= 20,
         AGEP <= 70) %>% 
  group_by(AGEP, sex, Major_Category, state) %>% 
  summarise(median = median(real_wage, na.rm = T),
            q90 = quantile(real_wage, 0.9, na.rm = T)) %>% 
  ungroup() %>% 
  mutate(model = "data")

bind_rows(ci_post_pred, data_intervals_median) %>% 
  mutate(plot_cat = paste0(str_sub(sex, 1, 1), " (", model, ")")) %>% 
  ggplot(., aes(x = AGEP)) +
  geom_line(aes(y = median, colour = plot_cat)) +
  geom_line(aes(y = q90, colour = plot_cat), linetype = "dashed") + 
  scale_colour_manual(values = c("#e31a1c", "#fd8d3c", "#08306b", "#2171b5"),
                      name = "group") +
  ggtitle("Posterior prediction (linear model) & real data: yearly wages (median & 90% quantile)") +
  xlab("age") +
  ylab("yearly wages in USD") +
  facet_grid(state ~ Major_Category)

The median values are roughly comparable with the ones in the real data, however our posterior sample shows higher variation than is found in real wage data. This is especially prominent in the education and health cohorts – it seems like our particular “hierarchical” implementation does not allow the Gaussian noise of the different major categories to vary a lot.

Now that we have posterior samples, we can use them to investigate some specific aspects of the wage data by deriving quantities from the posterior sample. One such derivation is to calculate the difference of predicted yearly wages for men and women in each panel.

ci_diff<- post_long %>% 
  spread(key = "sex", value = "predicted_earning") %>% 
  mutate(diff_mf = male - female) %>%
  group_by(AGEP, Major_Category, state) %>%
  summarise(q_lower = quantile(diff_mf, 0.025, names = F),
            median = quantile(diff_mf, 0.5, names = F),
            q_upper = quantile(diff_mf, 0.975, names = F))

ci_diff %>%
  ggplot(., aes(x = AGEP)) +
  geom_line(aes(y = median)) +
  geom_line(aes(y = q_lower), linetype = "dashed") +
  geom_line(aes(y = q_upper), linetype = "dashed") +
  ggtitle("Posterior prediction (linear model): wage difference (male - female; central 95% interval)") +
  xlab("age") +
  ylab("difference in yearly wages [USD]") +
  facet_grid(state ~ Major_Category)

It seems like there are disparities between the panels when looking at these wage differences. However, as we have alluded to earlier, this might be entirely due to the fact that the wage levels themselves also vary across the panels.

We can see if that is the case by investigating the wage ratio instead of the difference. This gives us a relative measure instead of an absolute one. A ratio like this is best visualised on the log2 scale:

ci_ratio <- post_long %>% 
  spread(key = "sex", value = "predicted_earning") %>% 
  mutate(ratio_mf = male / female) %>% 
  group_by(AGEP, Major_Category, state) %>% 
  summarise(q_lower = quantile(ratio_mf, 0.025, names = F),
            median = quantile(ratio_mf, 0.5, names = F),
            q_upper = quantile(ratio_mf, 0.975, names = F)) 

ci_ratio %>% 
  ggplot(., aes(x = AGEP)) +
  geom_line(aes(y = log2(median))) +
  geom_line(aes(y = log2(q_lower)), linetype = "dashed") +
  geom_line(aes(y = log2(q_upper)), linetype = "dashed") +
  ggtitle("Posterior prediction (linear model): log2 ratio of wages (male / female; central 95% interval)") +
  xlab("age") +
  ylab("log2 ratio of yearly wages") +
  facet_grid(state ~ Major_Category) 

Unfortunately (in terms of interesting outcomes), it does not seem like we can identify any group differences in the relative gender wage gap trend with this linear model: the slope and location of the time series are practically identical across all six panels.

As unexciting as it may be, this happenstance makes it easier to provide a summary of the predicted wage ratios for workers of different ages by simply taking the mean across all six panels:

ci_ratio %>% 
  filter(AGEP %in% seq(26, 50, 3)) %>% 
  group_by(AGEP) %>% 
  summarise(mean_lower = mean(q_lower),
            mean_50 = mean(median),
            mean_upper = mean(q_upper)) %>% 
  mutate(mean_lower = paste0("female ", round(1 / mean_lower, 2), "x higher"),
         mean_50 = paste0("male ", round(mean_50, 2), "x higher"),
         mean_upper = paste0("male ", round(mean_upper, 2), "x higher")) %>% 
  knitr::kable(., col.names = c("age", "2.5%", "50%", "97.5%"), align = "rccc")

age	2.5%	50%	97.5%
26	female 10.4x higher	male 1.26x higher	male 16.55x higher
29	female 9.76x higher	male 1.33x higher	male 17.57x higher
32	female 9.85x higher	male 1.39x higher	male 18.1x higher
35	female 9.3x higher	male 1.45x higher	male 19.26x higher
38	female 8.96x higher	male 1.51x higher	male 19.4x higher
41	female 8.62x higher	male 1.55x higher	male 21.36x higher
44	female 8.33x higher	male 1.63x higher	male 21.28x higher
47	female 7.82x higher	male 1.76x higher	male 22.14x higher
50	female 7.6x higher	male 1.79x higher	male 22.91x higher

This table shows the trends shown from the previous figure on a more tangible scale than the logarithmic scale: as men’s relative wages grow larger with age, the credibility intervals of wage ratios shift accordingly.